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ABSTRACT 

A  large  class  of  problems  in  engineering  mechanics  involves  a  so-called  “complementarity” 
relationship  representing  the  orthogonality  of  two  sign-constrained  vectors.  Typical  instances 
are  plasticity  laws  and  contact-like  conditions.  Eor  state  problems,  the  formulation  leads  to  a 
mixed  complementarity  problem  (MCP)  whereas  in  synthesis  (e.g.  minimum  weight  design) 
or  identification  problems,  a  mathematical  program  with  equilibrium  constraints  (MPEC)  is 
obtained.  The  aim  of  this  paper  is  two-fold.  Firstly,  it  describes,  through  two  typical  models, 
how  some  important  engineering  mechanics  problems  can  be  formulated  elegantly  and 
naturally  as  either  an  MCP  or  an  MPEC.  Secondly,  it  describes  a  powerful  computer-oriented 
environment  for  constructing  and  solving  these  mathematical  programming  problems,  with 
features  such  as  sparsity  and  automatic  differentiation  facilities  being  transparently  accessible. 
This  involves  the  use  of  the  modeling  language  GAMS  (an  acronym  for  General  Algebraic 
Modeling  System)  and  its  associated  mathematical  programming  solvers  (e.g.  the  industry 
standard  MCP  solver  PATH).  A  simple  generic  model  suitable  for  solving  the  state  problem 
for  trusses  is  used  to  clarify  the  syntax  of  GAMS  models  and  to  illustrate  the  ease  with  which 
they  can  be  built  and  solved. 
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INTRODUCTION 


Complementarity,  namely  the  requirement  that  two  sign-constrained  vectors  are  orthogonal,  is 
a  typical  and  recurrent  mathematical  structure  of  many  state,  design  and  inverse  problems  in 
nonlinear  mechanics.  This  was  first  recognized  in  the  late  1960s  by  Maier  whose  seminal 
work  led  some  years  later  to  the  NATO  conference  “Engineering  Plasticity  by  Mathematical 
Programming”  (Cohn  and  Maier  1979),  with  important  contributions  by  prominent 
researchers  from  both  engineering  and  mathematical  programming  communities.  The 
proceedings  of  that  workshop  still  represent  a  valuable  source  of  information  on  the  elegant 
and  powerful  framework  provided  by  mathematical  programming,  in  particular 
complementarity,  to  discrete  plasticity.  More  recent  sources  of  reference,  specifically  on  some 
engineering  and  economic  applications  of  complementarity  problems,  are  the  review  paper  by 
Ferris  and  Pang  (1997a)  and  the  proceedings  of  the  first  “International  Conference  on 
Complementarity  Problems”  (Ferris  and  Pang  1997b). 

To  engineers,  the  study  and  application  of  complementarity  notions  in  mechanics  should  have 
a  twofold  appeal:  a  refined  mathematical  formalism  rich  in  useful  theoretical  results  and  a 
wealth  of  efficient  and  robust  numerical  algorithms.  Unfortunately,  the  application  of  such 
concepts  has  been  sporadic  and  below  expectation.  For  instance,  plasticity  problems 
(involving  a  set  of  complementarity  conditions  between  yield  functions  and  plastic  strains)  are 
still  largely  solved  through  the  iterative  use  of  linear  solvers  when  a  complementarity 
formulation  would  automatically  choose  which  inequalities  to  satisfy  as  equations. 

The  motivation  of  this  paper  is  to  show  how  modeling  systems  and  their  associated 
mathematical  programming  solvers  can  help  with  model  building  and  solution  of  a  number  of 
important  engineering  mechanics  problems,  all  characterized  by  the  presence  of 
complementarity  conditions.  In  particular,  we  adopt  the  well-known  governing  equations  of  a 
simple  holonomic  (path-independent  or  nonlinear  elastic)  elastoplastic  discrete  model  (a)  to 
show  how  a  state  problem  can  be  formulated  as  a  mixed  complementarity  problem  (MCP);  (b) 
to  formulate  a  minimum  weight  problem  as  an  instance  of  a  mathematical  program  with 
equilibrium  (complementarity  in  our  case)  constraints  (MPEG);  and  (c)  to  describe  how  the 
modeling  system  GAMS,  an  acronym  for  General  Algebraic  Modeling  System  (Brooke  et  al. 
1992),  can  be  used  to  model  and  solve,  using  the  industry  standard  MCP  solver  PATH  (Dirkse 
and  Ferris  1995a),  a  simple  example  of  the  state  problem  for  elastoplastic  trusses. 


GOVERNING  RELATIONS  FOR  DISCRETE  HOLONOMIC  PLASTICITY 

We  refer  to  a  suitably  space-discretized  strucmral  system,  the  constituents  (finite  elements)  of 
which  obey  holonomic  plasticity  laws.  The  governing  relations  for  the  whole  structure  can  be 
elegantly  expressed  through  generalized  variables  (see  e.g.  Cohn  and  Maier  1979)  as  follows: 


/  =  C^x,  q  =  Cv, 
q  =  e  +  p,  x  =  Se,  p  =  Nz, 
w  =  -N^x  +  Hz  +  r>0,  z  >  0,  w^z  =  0. 


(1,2) 

(3,4,5) 

(6) 


Vector  and  matrix  quantities  represent  the  unassembled  contributions  of  corresponding 
elemental  entities  as  concatenated  vectors  and  block-diagonal  matrices,  respectively.  For  a 
structure  with  d  degrees  of  freedom  and  m  member  generalized  quantities,  Eqn.  1  expresses 
equilibrium,  through  compatibility  matrix  Ce  ,  between  the  nodal  loads  /  G  SI''  and  the 
natural  stresses  a;  e  Si” .  Eqn.  2  describes  linear  compatibility  of  strains  g  e  SH”  with  the 
nodal  displacements  v  6  91*^ .  Relations  3-6  embody  the  holonomic  constitutive  laws.  The 
additivity  of  elastic  ee  91"  and  plastic  pe  9^"  strains  is  given  by  Eqn.  3.  Linear  elasticity  is 

represented  in  Eqn.  4,  where  5g91"""'  collects  unassembled  element  stiffnesses.  Plastic 
strains  p  are  defined  in  Eqn.  5  by  an  associated  flow  rule  in  term  of  the  plastic  multipliers 
z  G  91^  (y  =  number  of  yield  functions)  through  the  matrix  of  outward  normals  N  G  9i"”‘^  to 
the  yield  surface.  Finally,  we  define  in  Eqn.  6  a  linear  (yield)  function  w(a:(z),  z)  :  9i^  — >  91’’ 
which  is  complementary  with  the  nonnegative  plastic  multiplier  vector  z  and  which 
accommodates,  through  H  G  9^’”'’’ ,  a  class  of  hardening  models  with  yield  limits  reSi’’ . 


COMPLEMENTARITY  MODELS 

Based  on  the  governing  relations  given  by  Eqns.  1-6,  we  are  now  in  a  position  to  formulate 
two  types  of  mathematical  problems  involving  complementarity.  As  representative  models, 
we  describe  a  standard  state  problem  cast  as  an  MCP  and  then  a  minimum  weight  synthesis 
problem  as  an  MPEC. 

State  Problem  (MCP) 

The  holonomic  state  problem  requires  the  calculation  of  the  state  variables  (a:,v,z)  for  a 
given  structure  (i.e.  for  specified  material  and  geometric  properties)  and  loading.  Since  we 
intend  to  use  the  GAMS  modeling  system  for  solving  this  problem,  we  adopt  a  “mixed” 
formulation  involving  both  static  (ac)  and  kinematic  variables  (v,  z)  (at  variance  with  the  usual 
approach  of  using  z  variables  only). 

After  some  obvious  substitutions,  the  problem  then  becomes  one  of  finding  (a:,  v,  z)  from  the 
following  relations: 


C"x-/  =  0, 

S’'a:-Cv  +  Vz  =  0, 

(7) 

w  s  -N^x  +  Hz  +  r>0,  z  >  0,  w^z  =  0, 

-  oo  <  (a:,  v)  <  +  oo. 

The  problem  given  by  Eqn.  7  is  an  example  of  a  general  MCP  (Dirkse  and  Ferris  1995a),  for 
which  it  is  required  to  find  z€3i"  for  given  lower  I  and  upper  bounds  u  (-o°  <  £  <  u  <  +  <>°) 
and  a  function  F :  91"  — >  91" ,  such  that  precisely  one  of  the  following  holds  for  each 

iG 


F^=0  and  i^<Zi<Ui, 


(8) 


F;>0  and  Zi=(i, 

Fi<0  and  Zi=u^. 

Our  MCP  (Eqn.  7)  can  be  solved  using  standard  methods  if  hardening  is  adopted.  In  certain 
instances,  however,  such  as  when  softening  laws  are  assumed,  there  is  no  guarantee  that  any 
of  known  algorithms  will  solve  the  problem  (see  e.g.  Tin-Loi  and  Ferris  1997). 


Minimum  Weight  Problem  (MPEC) 

The  minimum  weight  problem  we  use  as  an  example  of  an  MPEC  was  first  formulated  by 
Kaneko  and  Maier  (1981)  and  later  revisited  by  Ferris  and  Tin-Loi  (1999a)  with  a  view 
towards  more  efficient  and  robust  solution  schemes,  especially  for  large-scale  structures.  The 
problem  can  be  described  briefly  as  follows.  Under  the  assumptions  of  a  fixed  topology  and 
specified  loads,  we  wish  to  minimize  the  volume  of  the  structure  under  the  additional 
constraints  that  certain  or  all  displacements  and  plastic  deformations  are  kept  within 
prescribed  serviceability  limits.  The  yield  limits  r,  stiffnesses  S  and  hardening  parameters  H 
of  the  constituent  members  of  the  structure  are  all  regarded  as  unknown  but  assumed  to  be 
(continuous)  functions  of  the  cross-sectional  areas  of  all  n  elements. 

For  simplicity  of  exposition,  assume  a  truss-like  structure  of  n  members,  for  which  the 
unknown  element  areas  are  collected  in  vector  a  £  91"  and  the  known  element  lengths  in 
vector  /  €  91" .  Assuming  that  explicit  expressions  for  member  stiffnesses  S(a)  and  hardening 
matrices  H{a) ,  in  terms  of  a  are  available,  the  minimum  volume  (weight)  problem  can  then 
be  formally  stated  as  the  following  constrained  optimization  problem  in  (a,  x,  v,  z)  '■ 

min  l^a 

subject  to  C^x  -  /  =  0, 

X  -  S{a)Cv  +  Sia)Nz  =  0, 
w  =  -N^x  +  Hia)z  +  r(a)>0,  z>0,  w^z  =  0, 

W 

-  V  <  V  <  V, 

z  ^  z, 

Uf  <a<a„, 

Ta  =  0, 


where  ve  01*'  is  a  vector  of  nonnegative  deflections  limits;  z€  91’’  is  a  vector  of  prescribed 
upper  bounds  on  plastic  multipliers  used  to  model  the  limited  ductility  of  the  members; 
Ui  £  91"  and  £  91"  are,  respectively,  lower  and  upper  bounds  on  the  cross-sectional  areas; 


and  T  €  is  a  technological  matrix  imposing  t  constraints  (e.g.  identical  areas  for  groups 
of  members)  on  the  design  variables  (areas). 

The  optimization  problem  given  by  Eqn.  9  is  a  special  case  of  an  MPEC  (Luo  et  al.  1997)  in 
which  the  equilibrium  system  takes  the  form  of  a  complementarity  condition.  MPECs  are 
much  harder  to  solve  than  MCPs.  Whilst  an  extensive  theory  of  first  and  second  order 
optimality  conditions  for  MPECs  has  been  developed,  still  relatively  little  is  known  about  the 
numerical  solution  of  practical,  large-scale  MPECs  likely  to  arise  in  realistic  applications.  The 
most  prominent  feature  of  an  MPEC,  and  one  that  distinguishes  it  from  a  standard  nonlinear 
program,  is  the  presence  of  complementarity  constraints.  These  constraints  classify  this  class 
of  mathematical  programs  as  a  nonlinear  disjunctive  (or  piecewise)  program  and  therefore 
carries  with  it  a  “combinatorial  curse”.  Recent  work  by  Dirkse  and  Ferris  (1999)  describes 
several  new  tools  for  modeling  MPECs  that  are  built  around  the  introduction  of  an  MPEC 
model  type  into  the  GAMS  language,  ready  to  be  linked  to  newly  developed  solvers.  We, 
however,  have  had  considerable  success  in  modeling  and  solving  MPECs  for  a  variety  of 
engineering  mechanics  problems  (Ferris  and  Tin-Loi  1998,  1999a,  1999b)  as  a  series  of 
nonlinear  programming  problems  using  the  GAMS  environment  and  its  associated  nonlinear 
programming  solver  CONOPT  (Drud  1994). 


MODELING  WITH  GAMS 

GAMS  is  a  high-level  modeling  language  specially  designed  to  facilitate  the  construction, 
solution  and  maintenance  of  large  and  complex  mathematical  programming  models.  It  is  a 
high  level  declarative  language  for  formulating  small  to  very  large  mathematical 
programming  models  using  simple  and  concise  algebraic  statements  which  mirror  the  actual 
mathematical  constructs  involved.  A  GAMS  model  is  transparent  to  both  human  and 
computer,  is  easily  modified  and  moved  across  different  computing  platforms  from  notebooks 
to  mainframes,  and  is  independent  of  the  solution  algorithm  of  the  mathematical  programming 
solvers.  It  not  only  frees  the  model  builder  from  the  burdens  imposed  by  the  solution  phase 
but  also  takes  over  the  steps  required  for  generation  of  the  model.  In  addition  to  providing 
simplicity  and  compactness  of  model  construction,  it  possesses  important  capabilities  such  as 
an  internal  efficient  sparse  data  representation  and  automatic  differentiation. 

A  number  of  mathematical  programming  problems  types  can  be  solved  via  GAMS.  In 
addition  to  the  MCP  problem  type,  other  available  model  types  are  LP  (linear  programming), 
NLP  (nonlinear  programming),  MIP  (mixed  integer  programming),  RMIP  (relaxed  mixed 
integer  programming),  MINLP  (mixed  integer  nonlinear  programming),  RMINLP  (relaxed 
mixed  integer  nonlinear  progranuning)  and  CNS  (constrained  nonlinear  systems).  GAMS  is 
continually  evolving  and  adapted  as  new  algorithms  and  problem  classes  have  been  explored. 
We  refer  the  interested  reader  to  the  extensive  GAMS  library  of  models  (from  such  diverse 
areas  as  economics,  chemical  engineering,  trade,  etc.)  accessible  from  the  GAMS  website 
(http :  /  /  WWW .  gams .  com),  and  to  Dirkse  and  Ferris  (1995b)  for  GAMS  models  of  MCPs. 

In  order  to  illustrate  the  typical  structure  and  syntax  of  a  GAMS  model,  we  list  in  Fig.  1  a 
simple  model  (state. gms)  suitable  for  the  large-scale  holonomic  analysis  of  elastoplastic 


trusses.  The  GAMS  file  is  written  using  a  standard  text  editor  and  executed  through  a  “gams 
state”  command.  Readers  versed  in  GAMS  will  recognize  the  sets,  variables,  etc. 
declarations,  while  those  not  familiar  with  GAMS  will  appreciate  the  concise  yet  descriptive 
style  and  will  also  immediately  recognize  the  parallel  to  the  MCP  given  by  Eqn.  7.  We  have 
purposely  separated  the  model  proper  from  its  input  data  which  is  inserted  at  compile  time 
through  the  $  include  state.dat  statement.  Note  also  the  (optional)  matching  of  free 
variables  to  equations  in  the  model  statement,  allowing  GAMS  to  check  that  there  are  the 
same  number  of  free  variables  as  equations. 


sets 

d 

'No. 

of 

structure  dof' 

m 

'No, 

of 

members ' 

y 

'No. 

of 

yield  functions  per  member 

alias  (y,yy); 

parameters 
f  (d) 

C  (m,  d) 

S  (m) 

N{m,y) 

H{m,y,y) 
r  (m,y) 

variables 

x{m)  'Member  stresses' 

v{d)  'Displacements'; 

positive  variables 

z(m,y)  'Plastic  multipliers'; 

equations 
eql (d) 
eq2 (m) 
eq3  {m,y)  ; 

eql(d)  ..  sum(m,C(m,d) *x(m) ) -f (d)  =e=  0; 

eq2  (m)  ,.  (1/S{m)  )  *x{m) -sum(d,C(m,d)  *v(d)  ) +s\im(y,N(m,y)  *z  (m,y)  )  =e=  0; 

eq3(m,y)  ..  -N(m,y)  *x{m) +s\im{yy,H(m,y,yy)  *z  (m,yy)  ) +r  (m,y)  =g=  0; 

model  state  /eql.v,  eq2.x,  eq3.z/; 

■$ include  state.dat 

solve  state  using  mcp; 

display  v.l,x.l,z.l; 

Figure  1 :  A  simple  GAMS  MCP  state  model  state .  gms 


'Load  vector' 
'Compatibility  matrix' 
'Member  stiffness' 
'Normal  matrix' 
'Hardening  matrix' 
'Yield  limits'; 


As  a  simple  and  specific  example  of  data  input,  consider  the  academic  three-bar  truss  shown 
in  Fig.  2.  Relevant  data  state.dat  in  GAMS  format  are  also  indicated.  The  compatibility 
matrix  for  this  simple  example  is  explicitly  entered  (rather  than  generated  as  it  would  be  for  a 
large  problem).  A  simple  noninteracting  hardening  matrix  (with  nonzero  diagonal  entries)  is 
assumed.  On  executing  this  GAMS  model  using  the  default  solver  PATH  a  displacement 
vector  V  =  (1.518, 0.642)  and  indication  that  only  bar  1  yields  (in  tension)  will  be  obtained. 


600 


400 


sets  d  /  dl*d2  / 

m  /  ml*m3  / 

y  /  yl*y2  /; 


f 

("dl ") 

=  400; 

f("d2")  =  600; 

c 

( "ml" ,  ' 

dl")  = 

0.6; 

c 

("ml",  ’ 

d2")  = 

GO 

O 

c 

( "m2" ,  ' 

d2")  = 

1; 

c 

( "m3  " ,  ■ 

dl")  = 

-0,6; 

c 

{ "m3  " ,  ’ 

d2")  = 

00 

o 

s 

("ml ") 

=  400; 

S("m2")  =  500; 

s 

( "m3 " ) 

=  S("ml 

.")  ; 

N 

(m,  "yl" 

)  =  1; 

N(m, "y2")  =  -1 

H 

(m,y,y) 

=  0.125*S(m)  ; 

r 

(m,y)  = 

500; 

Figure  2;  Three-bar  truss  and  associated  state .  dat  input 


CONCLUSIONS 

A  large  number  of  important  problems  in  nonlinear  mechanics  (e.g.  plasticity  and  those  with 
contact-like  conditions)  involve  complementarity  relationships.  Indeed,  as  shown  in  this 
paper,  the  most  natural,  elegant  and  powerful  method  of  tackling  these  problems  is  often  as 
mathematical  programming  problems  involving  the  complementarity  relations  explicitly.  The 
two  main  problem  classes  (MCP  and  MPEC)  which  arise  can  both  be  modeled  and  solved 
within  the  GAMS  modeling  environment,  using  its  industry  standard  solvers. 

We  argue,  by  illustrating  how  easily  MCPs  and  MPECs  can  be  formulated  and  modeled  for  a 
specific  instance  of  holonomic  plasticity,  that  modeling  systems  can  provide  the  impetus 
required  for  wider  use  of  mathematical  programming  methods  in  solving  practical  problems 
involving  complementarity  conditions  in  engineering  mechanics.  Hopefully,  such  work  will 
also  lead  to  increased  synergetic  interaction  between  modelers  and  algorithm  developers. 
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